3. Full SEM of NCRMP data
Sara Devi Swaminathan
2023-08-22
#Set working directory to project directory and load packages
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
library(knitr)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.1.8
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(Hmisc)## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(brms)## Loading required package: Rcpp
## Loading 'brms' package (version 2.18.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
##
## The following object is masked from 'package:survival':
##
## kidney
##
## The following object is masked from 'package:stats':
##
## arlibrary(parallel)
library(ggridges)
library(loo)## This is loo version 2.5.1
## - Online documentation and vignettes at mc-stan.org/loo
## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.library(bayesplot)## This is bayesplot version 1.10.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme settingLoad data and set instructions for sampler
ncrmp <- read_csv("data/outputs/ncrmp_simple.csv")## New names:
## Rows: 1637 Columns: 49
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (5): SampleID, REGION, SUB_REGION_NAME, rugosityMethod, SEASON dbl (44): ...1,
## totalBiomass, coralAssociated, herb, otherFish, YEAR, MONTH,...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`logit <- function(mu){
logit.mu <- log(mu/(1-mu))
}
#-Sampler instructions-#
n_cores <- detectCores()
n_chains <- 2
n_iter <- 4000
n_warmup <- 2000Build individual models
SCTLD-susceptible coral
################################################################################
#-Build individual models-######################################################
################################################################################
#-Susceptible coral-############################################################
suscCoral_mod <- bf(propSusc ~ 1 +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
phi ~ 1 +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
zi ~ 1 +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
get_prior(suscCoral_mod, data = ncrmp)## prior class coef group resp dpar nlpar lb
## (flat) b
## (flat) b sMAX_DEPTH_1
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd SUB_REGION_NAME 0
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME 0
## student_t(3, 0, 2.5) sd YEAR 0
## student_t(3, 0, 2.5) sd Intercept YEAR 0
## student_t(3, 0, 2.5) sds 0
## student_t(3, 0, 2.5) sds s(MAX_DEPTH) 0
## (flat) b phi
## (flat) b sMAX_DEPTH_1 phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd SUB_REGION_NAME phi 0
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME phi 0
## student_t(3, 0, 2.5) sd YEAR phi 0
## student_t(3, 0, 2.5) sd Intercept YEAR phi 0
## student_t(3, 0, 2.5) sds phi 0
## student_t(3, 0, 2.5) sds s(MAX_DEPTH) phi 0
## (flat) b zi
## (flat) b sMAX_DEPTH_1 zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd SUB_REGION_NAME zi 0
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME zi 0
## student_t(3, 0, 2.5) sd YEAR zi 0
## student_t(3, 0, 2.5) sd Intercept YEAR zi 0
## student_t(3, 0, 2.5) sds zi 0
## student_t(3, 0, 2.5) sds s(MAX_DEPTH) zi 0
## ub source
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)propSusc_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(propSusc)) %>%
mutate(logitmean = logit(mean))
print(paste("The logit-transformed propSusc intercept is", round(propSusc_mean[1,3], 2)))## [1] "The logit-transformed propSusc intercept is -2.06"# "The logit-transformed propSusc intercept is -2.06"
susc_prior <- c(prior(normal(-2.1, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(normal(0, 2), class = "sds"),
prior(exponential(2), class = "sd"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(exponential(2), class = "sd", dpar = "zi"),
prior(normal(0, 2), class = "sds", dpar = "zi"))
susc_prior$resp <- "propSusc"SCTLD-resistant coral
#-Resistant coral-##############################################################
resistantCoral_mod <- bf(propResistant ~ 1 +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
phi ~ 1 +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
zi ~ 1 +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
get_prior(resistantCoral_mod, data = ncrmp) ## prior class coef group resp dpar nlpar lb
## (flat) b
## (flat) b sMAX_DEPTH_1
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd SUB_REGION_NAME 0
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME 0
## student_t(3, 0, 2.5) sd YEAR 0
## student_t(3, 0, 2.5) sd Intercept YEAR 0
## student_t(3, 0, 2.5) sds 0
## student_t(3, 0, 2.5) sds s(MAX_DEPTH) 0
## (flat) b phi
## (flat) b sMAX_DEPTH_1 phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd SUB_REGION_NAME phi 0
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME phi 0
## student_t(3, 0, 2.5) sd YEAR phi 0
## student_t(3, 0, 2.5) sd Intercept YEAR phi 0
## student_t(3, 0, 2.5) sds phi 0
## student_t(3, 0, 2.5) sds s(MAX_DEPTH) phi 0
## (flat) b zi
## (flat) b sMAX_DEPTH_1 zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd SUB_REGION_NAME zi 0
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME zi 0
## student_t(3, 0, 2.5) sd YEAR zi 0
## student_t(3, 0, 2.5) sd Intercept YEAR zi 0
## student_t(3, 0, 2.5) sds zi 0
## student_t(3, 0, 2.5) sds s(MAX_DEPTH) zi 0
## ub source
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)propResistant_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(propResistant)) %>%
mutate(logitmean = logit(mean))
print(paste("The logit-transformed propResistant intercept is", round(propResistant_mean[1,3], 2)))## [1] "The logit-transformed propResistant intercept is -3.69"# "The logit-transformed propResistant intercept is -3.69"
resistant_prior <- c(prior(normal(-3.7, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(normal(0, 2), class = "sds"),
prior(exponential(2), class = "sd"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(exponential(2), class = "sd", dpar = "zi"),
prior(normal(0, 2), class = "sds", dpar = "zi"))
resistant_prior$resp <- "propResistant"Cyanobacteria
#-Cyanobacteria-################################################################
cyano_mod <- bf(cyanobacteria ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
phi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
zi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
cyano_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(cyanobacteria)) %>%
mutate(logit_mean = logit(mean))
print(paste("The logit-transformed cyanobacteria intercept is", round(cyano_mean[1,3],2)))## [1] "The logit-transformed cyanobacteria intercept is -1.47"#""The logit-transformed cyanobacteria intercept is -1.47"
get_prior(cyano_mod, data = ncrmp) ## prior class coef group resp
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## logistic(0, 1) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## phi default
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi default
## phi 0 default
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 default
## phi 0 (vectorized)
## zi default
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi default
## zi 0 default
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 default
## zi 0 (vectorized)cyano_prior <- c(prior(normal(-1.5, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(exponential(2), class = "sd"),
prior(normal(0, 2), class = "sds"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(exponential(2), class = "sd", dpar = "zi"),
prior(normal(0, 2), class = "sds", dpar = "zi"))
cyano_prior$resp <- "cyanobacteria" CCA
#-CCA-##########################################################################
cca_mod <- bf(CCA ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
zi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
phi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
CCA_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(CCA)) %>%
mutate(logit_mean = logit(mean))
print(paste("The logit-transformed CCA intercept is", round(CCA_mean[1,3], 2)))## [1] "The logit-transformed CCA intercept is -2.99"# "The logit-transformed CCA intercept is -2.99"
get_prior(cca_mod, data = ncrmp)## prior class coef group resp
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## logistic(0, 1) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## phi default
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi default
## phi 0 default
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 default
## phi 0 (vectorized)
## zi default
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi default
## zi 0 default
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 default
## zi 0 (vectorized)cca_brm_priors <- c(prior(normal(-3, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(exponential(2), class = "sd"),
prior(normal(0, 2), class = "sds"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(exponential(2), class = "sd", dpar = "zi"),
prior(normal(0, 2), class = "sds", dpar = "zi"))
cca_prior <- cca_brm_priors
cca_prior$resp <- "CCA"Macroalgae
#-Macroalgae-###################################################################
macro_mod <- bf(macroalgae ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
zi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
phi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
macroalgae_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(macroalgae)) %>%
mutate(logit_mean = logit(mean))
print(paste("The logit-transformed macroalgae intercept is", round(macroalgae_mean[1,3], 2)))## [1] "The logit-transformed macroalgae intercept is -0.95"# "The logit-transformed macroalgae intercept is -0.95"
get_prior(macro_mod, data = ncrmp)## prior class coef group resp
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## logistic(0, 1) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## phi default
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi default
## phi 0 default
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 default
## phi 0 (vectorized)
## zi default
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi default
## zi 0 default
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 default
## zi 0 (vectorized)macro_prior <- c(prior(normal(-1, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(exponential(2), class = "sd"),
prior(normal(0, 2), class = "sds"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(exponential(2), class = "sd", dpar = "zi"),
prior(normal(0, 2), class = "sds", dpar = "zi"))
macro_prior$resp <- "macroalgae" Fire coral
#-Fire coral-###################################################################
fire_mod <- bf(fireCoral ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
zi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
phi ~ 1 + rugosity_GC + mean_SR_rugosity +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
fireCoral_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(fireCoral)) %>%
mutate(logit_mean = logit(mean))
print(paste("The logit-transformed fireCoral intercept is", round(fireCoral_mean[1,3], 2)))## [1] "The logit-transformed fireCoral intercept is -5.16"# "The logit-transformed fireCoral intercept is -5.16
get_prior(fire_mod, data = ncrmp)## prior class coef group resp
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## logistic(0, 1) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## phi default
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi (vectorized)
## phi default
## phi 0 default
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 (vectorized)
## phi 0 default
## phi 0 (vectorized)
## zi default
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi (vectorized)
## zi default
## zi 0 default
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 (vectorized)
## zi 0 default
## zi 0 (vectorized)fire_prior <- c(prior(normal(-5.2, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(normal(0, 2), class = "sds"),
prior(exponential(2), class = "sd"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(exponential(2), class = "sd", dpar = "zi"),
prior(normal(0, 2), class = "sds", dpar = "zi"))
fire_prior$resp <- "fireCoral"Rugosity
#-Rugosity-#####################################################################
rug_mod <- bf(logNewRug ~ 1 + susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
rugosityMethod + (1|SUB_REGION_NAME) + (1|YEAR),
sigma ~ 1 + susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
rugosityMethod + (1|SUB_REGION_NAME) + (1|YEAR),family = student())
get_prior(rug_mod, data = ncrmp) ## prior class coef group resp
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosityMethodWTD_RUG
## (flat) b susc_GC
## student_t(3, -1, 2.5) Intercept
## gamma(2, 0.1) nu
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## (flat) b
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosityMethodWTD_RUG
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 1 default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## sigma default
## sigma (vectorized)
## sigma (vectorized)
## sigma (vectorized)
## sigma (vectorized)
## sigma (vectorized)
## sigma default
## sigma 0 default
## sigma 0 (vectorized)
## sigma 0 (vectorized)
## sigma 0 (vectorized)
## sigma 0 (vectorized)logNewRug_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(logNewRug))
print(paste("The logNewRug intercept is", round(logNewRug_mean[1,2], 2)))## [1] "The logNewRug intercept is -0.95"# "The logNewRug intercept is -0.95
rug_prior <- c(prior(normal(-1, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(exponential(1), class = "sd"),
prior(normal(-1, 1), class = "Intercept", dpar = "sigma"),
prior(exponential(1), class = "sd", dpar = "sigma"))
rug_prior$resp <- "logNewRug"Coral-associated fish
#-Coral associated fish-########################################################
coralAssoc_mod<- bf(coralAssociated ~ 1 +
rugosity_GC+ mean_SR_rugosity + rugosityMethod +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
mean_SR_cyano + cyano_GC +
mean_SR_cca + cca_GC +
mean_SR_macro + macro_GC +
mean_SR_fire + fire_GC +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
shape ~ 1 +
rugosity_GC+ mean_SR_rugosity + rugosityMethod +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
mean_SR_cyano + cyano_GC +
mean_SR_cca + cca_GC +
mean_SR_macro + macro_GC +
mean_SR_fire + fire_GC +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
hu ~ 1 +
rugosity_GC+ mean_SR_rugosity + rugosityMethod +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
mean_SR_cyano + cyano_GC +
mean_SR_cca + cca_GC +
mean_SR_macro + macro_GC +
mean_SR_fire + fire_GC +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = hurdle_gamma(link = "log", link_shape = "log", link_hu = "logit"))
get_prior(coralAssoc_mod, data = ncrmp)## prior class coef group resp
## (flat) b
## (flat) b cca_GC
## (flat) b cyano_GC
## (flat) b fire_GC
## (flat) b macro_GC
## (flat) b mean_SR_cca
## (flat) b mean_SR_cyano
## (flat) b mean_SR_fire
## (flat) b mean_SR_macro
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b rugosityMethodWTD_RUG
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0.4, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b cca_GC
## (flat) b cyano_GC
## (flat) b fire_GC
## (flat) b macro_GC
## (flat) b mean_SR_cca
## (flat) b mean_SR_cyano
## (flat) b mean_SR_fire
## (flat) b mean_SR_macro
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b rugosityMethodWTD_RUG
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## logistic(0, 1) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b cca_GC
## (flat) b cyano_GC
## (flat) b fire_GC
## (flat) b macro_GC
## (flat) b mean_SR_cca
## (flat) b mean_SR_cyano
## (flat) b mean_SR_fire
## (flat) b mean_SR_macro
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b rugosityMethodWTD_RUG
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## hu default
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu (vectorized)
## hu default
## hu 0 default
## hu 0 (vectorized)
## hu 0 (vectorized)
## hu 0 (vectorized)
## hu 0 (vectorized)
## hu 0 default
## hu 0 (vectorized)
## shape default
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape default
## shape 0 default
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 default
## shape 0 (vectorized)coralAssociated_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(coralAssociated)) %>%
mutate(logmean = log(mean))
print(paste("The log-transformed coralAssociated intercept is", round(coralAssociated_mean[1,3], 2)))## [1] "The log-transformed coralAssociated intercept is 2.31"# "The log-transformed coralAssociated intercept is 2.31"
coralAssoc_prior <- c(prior(normal(2.3, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(normal(0, 2), class = "sds"),
prior(exponential(1), class = "sd"),
prior(normal(2.3, 1), class = "Intercept", dpar = "shape"),
prior(exponential(1), class = "sd", dpar = "shape"),
prior(normal(0, 1), class = "b", dpar = "shape"))
coralAssoc_prior$resp <- "coralAssociated"Herbivorous fish
#-Herbivorous fishes-###########################################################
herb_mod <- bf(herb ~ 1 + mean_SR_rugosity + rugosity_GC +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
mean_SR_cyano + cyano_GC +
mean_SR_cca + cca_GC +
mean_SR_macro + macro_GC +
mean_SR_fire + fire_GC +
(1 | SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
shape ~ 1 + mean_SR_rugosity + rugosity_GC +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
mean_SR_cyano + cyano_GC +
mean_SR_cca + cca_GC +
mean_SR_macro + macro_GC +
mean_SR_fire + fire_GC +
(1 | SUB_REGION_NAME) + (1|YEAR)+ s(MAX_DEPTH),
family = hurdle_gamma(link = "log"))
get_prior(herb_mod, data = ncrmp)## prior class coef group resp
## (flat) b
## (flat) b cca_GC
## (flat) b cyano_GC
## (flat) b fire_GC
## (flat) b macro_GC
## (flat) b mean_SR_cca
## (flat) b mean_SR_cyano
## (flat) b mean_SR_fire
## (flat) b mean_SR_macro
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## beta(1, 1) hu
## student_t(3, -0.6, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b cca_GC
## (flat) b cyano_GC
## (flat) b fire_GC
## (flat) b macro_GC
## (flat) b mean_SR_cca
## (flat) b mean_SR_cyano
## (flat) b mean_SR_fire
## (flat) b mean_SR_macro
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## 0 1 default
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## shape default
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape default
## shape 0 default
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 default
## shape 0 (vectorized)herb_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(herb)) %>%
mutate(logmean = log(mean))
print(paste("The log-transformed herb intercept is", round(herb_mean[1,3], 2)))## [1] "The log-transformed herb intercept is -0.39"# "The log-transformed herb intercept is -0.39"
herb_prior <- c(prior(normal(-0.4, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(normal(0, 2), class = "sds"),
prior(exponential(1), class = "sd"),
prior(beta(1, 1), class = "hu"),
prior(normal(-0.4, 1), class = "Intercept", dpar = "shape"),
prior(exponential(1), class = "sd", dpar = "shape"),
prior(normal(0, 1), class = "b", dpar = "shape"))
herb_prior$resp <- "herb"Other fish
#-Other fish-###################################################################
otherFish_mod <- bf(otherFish ~ 1 +
rugosity_GC+ mean_SR_rugosity + rugosityMethod +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
mean_SR_cyano + cyano_GC +
mean_SR_cca + cca_GC +
mean_SR_macro + macro_GC +
mean_SR_fire + fire_GC +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
shape ~ 1 +
rugosity_GC+ mean_SR_rugosity + rugosityMethod +
susc_GC + mean_SR_susc_cover +
resistant_GC + mean_SR_resistant_cover +
mean_SR_cyano + cyano_GC +
mean_SR_cca + cca_GC +
mean_SR_macro + macro_GC +
mean_SR_fire + fire_GC +
(1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
family = Gamma(link = "log"))
get_prior(otherFish_mod, data = ncrmp)## prior class coef group resp
## (flat) b
## (flat) b cca_GC
## (flat) b cyano_GC
## (flat) b fire_GC
## (flat) b macro_GC
## (flat) b mean_SR_cca
## (flat) b mean_SR_cyano
## (flat) b mean_SR_fire
## (flat) b mean_SR_macro
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b rugosityMethodWTD_RUG
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 1, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## (flat) b
## (flat) b cca_GC
## (flat) b cyano_GC
## (flat) b fire_GC
## (flat) b macro_GC
## (flat) b mean_SR_cca
## (flat) b mean_SR_cyano
## (flat) b mean_SR_fire
## (flat) b mean_SR_macro
## (flat) b mean_SR_resistant_cover
## (flat) b mean_SR_rugosity
## (flat) b mean_SR_susc_cover
## (flat) b resistant_GC
## (flat) b rugosity_GC
## (flat) b rugosityMethodWTD_RUG
## (flat) b sMAX_DEPTH_1
## (flat) b susc_GC
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd SUB_REGION_NAME
## student_t(3, 0, 2.5) sd Intercept SUB_REGION_NAME
## student_t(3, 0, 2.5) sd YEAR
## student_t(3, 0, 2.5) sd Intercept YEAR
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(MAX_DEPTH)
## dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## shape default
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape (vectorized)
## shape default
## shape 0 default
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 (vectorized)
## shape 0 default
## shape 0 (vectorized)otherFish_mean <- ncrmp %>%
dplyr::group_by(SUB_REGION_NAME) %>%
dplyr::summarize(mean = mean(otherFish)) %>%
mutate(logmean = log(mean))
print(paste("The log-transformed otherFish intercept is", round(otherFish_mean[1,3], 2)))## [1] "The log-transformed otherFish intercept is 3.49"# "The log-transformed otherFish intercept is 3.49
otherFish_prior <- c(prior(normal(3.5, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(normal(0, 2), class = "sds"),
prior(exponential(1), class = "sd"),
prior(normal(3.1, 1), class = "Intercept", dpar ="shape"),
prior(normal(0, 1), class = "b", dpar = "shape"),
prior(exponential(1), class = "sd", dpar = "shape"))
otherFish_prior$resp <- "otherFish"Build and run full SEM
################################################################################
#-Build full SEM-###############################################################
################################################################################
full_priors <- c( susc_prior,
resistant_prior,
cyano_prior,
cca_prior,
macro_prior,
fire_prior,
rug_prior,
coralAssoc_prior,
herb_prior,
otherFish_prior)
full_priors## prior class coef group resp dpar nlpar lb ub
## normal(-2.1, 1) Intercept propSusc <NA> <NA>
## normal(0, 1) b propSusc <NA> <NA>
## normal(0, 2) sds propSusc <NA> <NA>
## exponential(2) sd propSusc <NA> <NA>
## logistic(0, 1) Intercept propSusc zi <NA> <NA>
## exponential(2) sd propSusc zi <NA> <NA>
## normal(0, 2) sds propSusc zi <NA> <NA>
## normal(-3.7, 1) Intercept propResistant <NA> <NA>
## normal(0, 1) b propResistant <NA> <NA>
## normal(0, 2) sds propResistant <NA> <NA>
## exponential(2) sd propResistant <NA> <NA>
## logistic(0, 1) Intercept propResistant zi <NA> <NA>
## exponential(2) sd propResistant zi <NA> <NA>
## normal(0, 2) sds propResistant zi <NA> <NA>
## normal(-1.5, 1) Intercept cyanobacteria <NA> <NA>
## normal(0, 1) b cyanobacteria <NA> <NA>
## exponential(2) sd cyanobacteria <NA> <NA>
## normal(0, 2) sds cyanobacteria <NA> <NA>
## logistic(0, 1) Intercept cyanobacteria zi <NA> <NA>
## exponential(2) sd cyanobacteria zi <NA> <NA>
## normal(0, 2) sds cyanobacteria zi <NA> <NA>
## normal(-3, 1) Intercept CCA <NA> <NA>
## normal(0, 1) b CCA <NA> <NA>
## exponential(2) sd CCA <NA> <NA>
## normal(0, 2) sds CCA <NA> <NA>
## logistic(0, 1) Intercept CCA zi <NA> <NA>
## exponential(2) sd CCA zi <NA> <NA>
## normal(0, 2) sds CCA zi <NA> <NA>
## normal(-1, 1) Intercept macroalgae <NA> <NA>
## normal(0, 1) b macroalgae <NA> <NA>
## exponential(2) sd macroalgae <NA> <NA>
## normal(0, 2) sds macroalgae <NA> <NA>
## logistic(0, 1) Intercept macroalgae zi <NA> <NA>
## exponential(2) sd macroalgae zi <NA> <NA>
## normal(0, 2) sds macroalgae zi <NA> <NA>
## normal(-5.2, 1) Intercept fireCoral <NA> <NA>
## normal(0, 1) b fireCoral <NA> <NA>
## normal(0, 2) sds fireCoral <NA> <NA>
## exponential(2) sd fireCoral <NA> <NA>
## logistic(0, 1) Intercept fireCoral zi <NA> <NA>
## exponential(2) sd fireCoral zi <NA> <NA>
## normal(0, 2) sds fireCoral zi <NA> <NA>
## normal(-1, 1) Intercept logNewRug <NA> <NA>
## normal(0, 1) b logNewRug <NA> <NA>
## exponential(1) sd logNewRug <NA> <NA>
## normal(-1, 1) Intercept logNewRug sigma <NA> <NA>
## exponential(1) sd logNewRug sigma <NA> <NA>
## normal(2.3, 1) Intercept coralAssociated <NA> <NA>
## normal(0, 1) b coralAssociated <NA> <NA>
## normal(0, 2) sds coralAssociated <NA> <NA>
## exponential(1) sd coralAssociated <NA> <NA>
## normal(2.3, 1) Intercept coralAssociated shape <NA> <NA>
## exponential(1) sd coralAssociated shape <NA> <NA>
## normal(0, 1) b coralAssociated shape <NA> <NA>
## normal(-0.4, 1) Intercept herb <NA> <NA>
## normal(0, 1) b herb <NA> <NA>
## normal(0, 2) sds herb <NA> <NA>
## exponential(1) sd herb <NA> <NA>
## beta(1, 1) hu herb <NA> <NA>
## normal(-0.4, 1) Intercept herb shape <NA> <NA>
## exponential(1) sd herb shape <NA> <NA>
## normal(0, 1) b herb shape <NA> <NA>
## normal(3.5, 1) Intercept otherFish <NA> <NA>
## normal(0, 1) b otherFish <NA> <NA>
## normal(0, 2) sds otherFish <NA> <NA>
## exponential(1) sd otherFish <NA> <NA>
## normal(3.1, 1) Intercept otherFish shape <NA> <NA>
## normal(0, 1) b otherFish shape <NA> <NA>
## exponential(1) sd otherFish shape <NA> <NA>
## source
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user
## user### UNCOMMENT BELOW TO RUN MODEL
# full_SEM <- brm(cyano_mod +
# rug_mod +
# cca_mod +
# fire_mod +
# macro_mod +
# otherFish_mod +
# coralAssoc_mod +
# herb_mod +
# resistantCoral_mod +
# suscCoral_mod +
# set_rescor(rescor = FALSE),
# prior = full_priors,
# data = ncrmp,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = "0",
# control = list(adapt_delta = 0.99))
## save to file
full_SEM <- readRDS("data/outputs/full_SEM.rds")
#saveRDS(full_SEM, "data/outputs/full_SEM.rds")Posterior predictive checks
SCTLD-susceptible coral
pp_check(full_SEM, ndraws = 100,resp = 'propSusc', type = "dens_overlay")+
labs(title = "Sctld-susceptible coral")
pp_check(full_SEM, ndraws = 100,resp = 'propSusc', type = "hist")+
labs(title = "Sctld-susceptible coral")
pp_check(full_SEM, ndraws = 100,resp = 'propSusc', type = "boxplot")+
labs(title = "Sctld-susceptible coral")
y <- t(posterior_predict(full_SEM, resp = 'propSusc')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$propSusc*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
SCTLD-resistant coral
pp_check(full_SEM, ndraws = 100,resp = 'propResistant', type = "dens_overlay")+
labs(title = "Sctld-resistant coral")
pp_check(full_SEM, ndraws = 100,resp = 'propResistant', type = "hist")+
labs(title = "Sctld-resistant coral")
pp_check(full_SEM, ndraws = 100,resp = 'propResistant', type = "boxplot")+
labs(title = "Sctld-resistant coral")
y <- t(posterior_predict(full_SEM, resp = 'propResistant')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$propResistant*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
cyanobacteria
pp_check(full_SEM, ndraws = 100,resp = 'cyanobacteria', type = "dens_overlay")+
labs(title = "cyanobacteria")
pp_check(full_SEM, ndraws = 100,resp = 'cyanobacteria', type = "hist")+
labs(title = "cyanobacteria")
pp_check(full_SEM, ndraws = 100,resp = 'cyanobacteria', type = "boxplot")+
labs(title = "cyanobacteria")
y <- t(posterior_predict(full_SEM, resp = 'cyanobacteria')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$cyanobacteria*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
CCA
pp_check(full_SEM, ndraws = 100,resp = 'CCA', type = "dens_overlay")+
labs(title = "cca")
pp_check(full_SEM, ndraws = 100,resp = 'CCA', type = "hist")+
labs(title = "cca")
pp_check(full_SEM, ndraws = 100,resp = 'CCA', type = "boxplot")+
labs(title = "cca")
y <- t(posterior_predict(full_SEM, resp = 'CCA')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$CCA*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
macroalgae
pp_check(full_SEM, ndraws = 100,resp = 'macroalgae', type = "dens_overlay")+
labs(title = "macroalgae")
pp_check(full_SEM, ndraws = 100,resp = 'macroalgae', type = "hist")+
labs(title = "macroalgae")
pp_check(full_SEM, ndraws = 100,resp = 'macroalgae', type = "boxplot")+
labs(title = "macroalgae")
y <- t(posterior_predict(full_SEM, resp = 'macroalgae')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$macroalgae*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
fire coral
pp_check(full_SEM, ndraws = 100,resp = 'fireCoral', type = "dens_overlay")+
labs(title = "fire coral")
pp_check(full_SEM, ndraws = 100,resp = 'fireCoral', type = "hist")+
labs(title = "fire coral")
pp_check(full_SEM, ndraws = 100,resp = 'fireCoral', type = "boxplot")+
labs(title = "fire coral")
y <- t(posterior_predict(full_SEM, resp = 'fireCoral')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$fireCoral*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
rugosity
pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "dens_overlay")+
labs(title = "rugosity")
pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "hist")+
labs(title = "rugosity")
pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "boxplot")+
labs(title = "rugosity")
pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "loo_pit_overlay")+
labs(title = "rugosity")## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
coral-associated fish
pp_check(full_SEM, ndraws = 100,resp = 'coralAssociated', type = "dens_overlay")+
labs(title = "coral-associated fish")
pp_check(full_SEM, ndraws = 100,resp = 'coralAssociated', type = "hist")+
labs(title = "coral-associated fish")
pp_check(full_SEM, ndraws = 100,resp = 'coralAssociated', type = "boxplot")+
labs(title = "coral-associated fish")
y <- t(posterior_predict(full_SEM, resp = 'coralAssociated')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$fireCoral*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
herbivorous fish
pp_check(full_SEM, ndraws = 100,resp = 'herb', type = "dens_overlay")+
labs(title = "herbivorous fish")
pp_check(full_SEM, ndraws = 100,resp = 'herb', type = "hist")+
labs(title = "herbivorous fish")
pp_check(full_SEM, ndraws = 100,resp = 'herb', type = "boxplot")+
labs(title = "herbivorous fish")
y <- t(posterior_predict(full_SEM, resp = 'herb')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$fireCoral*100)
length(y_int) ## [1] 1637ncol(yrep_int)## [1] 1637ppc_bars(y_int, yrep_int)
ppc_rootogram(y_int, yrep_int)
other fish
pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "dens_overlay")+
labs(title = "otherFish")
pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "hist")+
labs(title = "otherFish")
pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "boxplot")+
labs(title = "otherFish")
pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "loo_pit_overlay")+
labs(title = "otherFish")## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.